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ABSTRACT 

Simulations predict that the dark matter halos of galaxies should have central cusps, while those 
inferred from observed galaxies do not have cusps. We demonstrate, using both linear perturbation theory 
and n-body simulations, that a disk bar, which should be ubiquitous in forming galaxies, can produce 
cores in cuspy CDM dark matter profiles within five bar orbital times. Simulations of forming galaxies 
suggest that one of Milky Way size could have a 10 kpc primordial bar; this bar will remove the cusp out 
to '^5 kpc in ~1.5 gigayears, while the disk only loses ~8% of its original angular momentum. An inner 
Lindblad-like resonance couples the rotating bar to orbits at all radii through the cusp, transferring the 
bar pattern angular momentum to the dark matter cusp, rapidly flattening it. This resonance disappears 
for profiles with cores and is responsible for a qualitative difference in bar driven halo evolution with and 
without a cusp. This bar induced evolution will have a profound effect on the structure and evolution of 
almost all galaxies. Hence, both to understand galaxy formation and evolution and to make predictions 
from theory it is necessary to resolve these dynamical processes. Unfortunately, correctly resolving these 
important dynamical processes in ab initio calculations of galaxy formation is a daunting task, requiring 
at least 4,000,000 halo particles using our SCF code, and probably requiring many times more particles 
when using noisier tree, direct summation, or grid based techniques, the usual methods employed in such 
calculations. 

Subject headings: galaxies: evolution — galaxies: halos — galaxies: kinematics and dynamics — 
cosmology: theory — dark matter 



1. INTRODUCTION 

The cold dark matter (CDM) model for structure for- 
mation has had a wide range of successes in explaining the 
observed universe. However, particularly at small scales, 
several vexing problems remain. Perhaps the most widely 
discussed shortcoming concerns the central density profiles 
of galaxies. Cold dark matter structure formation simula- 
tions predict a universal cuspy halo profile (Navarro et al. 
1997, hereafter NEW). This profile, p cx r-T(l + r/r,)^-^ 
or p oc r~'^ {l+{r /rs)^~'^)~^ with 7 = 1, was first presented 
by NEW based on a suite of coUisionless n-body simula- 
tions with different initial density fluctuation spectra and 
cosmological parameters. More recent work debates the 
value of 7 (Moore et al. 1998; Jing and Suto 2000), but 
most estimates are in the range 1 < 7 < 1.5. Even before 
the discovery of a universal density profile, several authors 
pointed out the apparent discrepancy between the cuspy 
central density profiles predicted by dissipationless numer- 
ical simulations and those inferred from the rotation curves 
of galaxies (Moore 1994; Flores and Primack 1994; Burk- 
ert 1995). Although some recent observational evidence 
claims a marginal consistency with such cuspy dark mat- 
ter profiles (van den Bosch and Swaters 2001), most do 
not (Cote et al. 2000; de Blok et al. 2001; Blais-Ouellette 
et al 2001). 

This led to a flurry of papers trying to explain the dis- 
crepancy. Given the numerous other successes of the CDM 
model most of these explanations involved altering the 
model only at small scales, which would only affect the 
cuspiness of the dark matter density profile. Some changed 
the fundamental nature of the dark matter particle itself: 



collisional dark matter (Spergel and Steinhardt 2000), de- 
caying dark matter (Cen 2001), fluid dark matter (Peebles 
2000), repulsive dark matter (Goodman 2000), and an- 
nihilating dark matter (Kaplinghat et al. 2000). Others 
changed the temperature of the dark particle from cold 
to warm to reduce the small scale power, i.e. warm dark 
matter (Hogan and Dalcanton 2000; Colin et al. 2000; Eke 
et al. 2001; Avila-Reese et al. 2001; Bode et al. 2001), while 
others altered the shape of the spectrum at small scales by 
changing inflation (Kamionkowski and Liddle 2000). At 
the most extreme, some authors suggested that the only 
solution was to change the nature of gravity itself (de Blok 
and McGaugh 1998). 

The simulations that predicted the existence of the cen- 
tral dark matter cusp (Navarro et al. 1997; Moore et al. 
1998; Jing and Suto 2000) only included the dark matter 
component. By using rotation curve analyses and light 
profiles, the same techniques used to infer the above men- 
tioned discrepancy, we know that the central regions of 
most galaxies are dominated by a baryonic component 
(van Albada and Sancisi 1986). The dynamics of the bary- 
onic component can be quite different than that of dark 
matter owing to its dissipative nature: baryons can shock 
and cool. Furthermore, even a dissipationless, stellar bary- 
onic component can have a markedly different geometry 
than pure dark matter since it originally formed from a 
dissipative component. These different dynamical proper- 
ties and the subsequent interactions of the baryonic com- 
ponent with the dark matter could greatly affect the dark 
matter structure. 

Since it is commonly assumed that the addition of a 
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baryonic component would cause the dark matter to adia- 
batically contract (Blumenthal et al. 1986) thereby exacer- 
bating the observed discrepancy, the compUcations caused 
by the baryonic component are usually dismissed. Even 
so, both Binncy ct al. (2001) and El-Zant ct al. (2001) dis- 
cuss mechanisms by which the baryonic component might 
produce a core. Much of the baryonic component will cool 
and likely form a disk (Fall and Efstathiou 1980). Disk 
bars and strong grand-design structure can strongly inter- 
act with the dark matter (Weinberg 1985: Hcrnriuist and 
Weinberg 1992) by transferring angular momentum to the 
spheroid-halo component. Quoting from Hernquist and 
Weinberg (1992) 

"Cosmological simulations of [dark] halo formation from 
plausible initial conditions generally yield halos whose den- 
sity structure resembles that of [a Hernquist profile] (Du- 
binski and Carlberg 1991). While somewhat uncertain, 
it would appear that this result is in conflict with ob- 
served rotation curves of at least some galaxies where the 
dark matter halos apparently have significant core regions 
with roughly constant density (e.g. Flores et al. (1993)). 
The simulations reported here [of rotating bars] provide, in 
principle, a mechanism for developing cores in halos long 
after they form." 

It is this mechanism that we explore in this paper. We 
present a brief description and simulations of the mech- 
anism in §2 and in §3 we discuss the dynamics in more 
detail along with other implications. 

2. THE MECHANISM 

In an equilibrium dark matter halo the average orbit is 
highly eccentric (e w 0.5 on average) and, therefore, has 
low angular momentum for its energy. Very little energy is 
required to remove the cusp; a minor addition of angular 
momentum provides enough circularization to exclude an 
eccentric orbit from the central region. The pattern of a 
rotating bar has enough angular momentum, if transferred 
to the cusp, to remove it altogether. Our basic picture is 
as follows: 

1. Start with a forming galaxy. CDM simulations 
predict that its dark matter distribution will be 
cuspy (e.g. an NEW profile). 

2. As the galaxy continues to form, smooth accretion 
of gas and stochastic merging and dissipation 

is followed by the settling of the baryons into a 
mostly gaseous, cold disk. 

3. As the disk becomes more massive and centrally 
concentrated, at some point it will overwhelm the 
halo support and form a bar as seen in simulations 
of galaxy formation that include a dissipative 
baryonic component (Katz and Gunn 1991; 
Steinmetz 1997). Since the disk is mostly gaseous 
and cold and the halo is quite disturbed, it makes 
it easy to form a very strong bar. Simulations 
suggest that for Milky Way sized galaxies early 
bars have semi-major axes approaching 10 kpc 
(Steinmetz 1997). 

4. The rotating bar pattern excites a gravitational 
wake in the dark matter. The wake trails the bar 



causing the bar to slow its rotation and transfer its 
angular momentum to the dark matter (Weinberg 
1985). This forms a core in the dark matter 
distribution. The 10 kpc primordial bar will 
remove the cusp out to -^5 kpc in ~1.5 gigayears. 

5. After the formation of the core early on, a classic 
stellar bar may form and experience low torque in 
the current epoch. 

It is easy to estimate that a strong bar can have impor- 
tant consequences for halo evolution. A non-axisymmetric 
bar force is dominated by its quadrupole component. A 
toy model for a rotating gravitational quadrupole is two 
masses in orbit at the same distance from the center of a 
galaxy but at opposite position angles. We can use the 
Chandrasekhar dynamical friction formula to estimate the 
time scale for transferring all of the bar's angular momen- 
tum to the halo. The answer is a few bar rotation peri- 
ods (sec Weinberg 1985). Since the total angular momen- 
tum in the bar approaches that of the dark matter halo 
within the optical radius, we conclude that the evolution of 
the dark matter halo induced by a bar can be significant. 
These estimates have been further refined by Hernquist 
and Weinberg (1992) and Debattista and Sellwood (1998, 
2000) using n-body simulations. 

The evolution of the halo may be estimated analyti- 
cally using a perturbation expansion of the coUisionlcss 
Boltzmann equation and a solution of the resulting ini- 
tial value problem. The zeroth-order solution specifies the 
equilibrium galax:y and the first-order solution determines 
the linear response of the galaxy to some perturbation. 
The second-order solution contains the first non-transient 
change in the underlying distribution. By taking the limit 
for the evolution time scale to be much larger than an 
orbital time, the transient contribution can be made arbi- 
trarily small. However, we are often interested in interme- 
diate time scales. Explicit comparisons with n-body sim- 
ulations shows that this approximation is acceptable even 
for a small number of orbital time scales. Mathematical 
details can be found in Weinberg (1985, 2001b). 

Figure la shows the perturbation theory prediction for 
the evolution of the density profile. The simulation is of 
a bar in an NFW profile with a concentration of 20. The 
disk contains half the mass within the bar radius and the 
bar contains 30% of the disk mass. We choose the coro- 
tation radius to be the NFW scale length and the bar 
radius is chosen to be 0.5 scale lengths. The results de- 
scribed below are only weakly sensitive to this choice. The 
bar figure is represented by a homogeneous ellipsoid with 
axis ratios of 1:0.5:0.05. We derive the bar force from the 
quadrupole term of the gravitational potential of the ellip- 
soid. Ignoring the higher order multipoles in the bar po- 
tential does not greatly change the results but does cause 
us to slightly underestimate the evolution (Hernquist and 
Weinberg 1992). 

The angular momentum transfer to the dark matter 
takes place at commensurabilities between the pattern 
speed and orbital frequencies and is dominated by a few 
strong resonances. The torque pushes the inner halo orbits 
to higher angular momentum and energy, removing gravi- 
tational support for the inner cusp. The halo expands and 
flattens the cusp. After a few bar rotation times (several 
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Fig. 1. — Evolution of an NFW profile with an embedded disk bar. Profiles are a time sequence labeled by the number of bar rotation 
times as predicted from linear perturbation theory (a) and from a n-body simulation (b). 



hundred million years) the central cusp is flattened, but 
not enough to match the observations. 

Although the pcrturbative approach yields more accu- 
rate results, free from the numerical noise inherent in n- 
body simulations, they are only valid in the linear regime. 
Hence, we perform n-body simulations using a parallel im- 
plementation with the Message Passing Interface (MPI), 
of the algorithm described in Weinberg (1999), a self con- 
sistent field (SCF) technique. This algorithm defines a set 
of orthogonal functions whose lowest-order member is the 
unperturbed profile itself. Each additional member in the 
series probes successively finer scales. Because all scales 
of interest here can be represented with a small number of 
degrees of freedom, the particle noise is low. 

The initial conditions are a Monte Carlo realization of 
the exact isotropic phase space distribution function for 
the NFW profile, determined by Eddington inversion (see 
Binney and Tremaine 1987, Chap. 4), using 4,000,000 
equal mass particles. As we discuss below this number 
of particles is sufficient to give a converged result. Fur- 
thermore, when simulated without a bar disturbance this 
realization does not evolve and remains in equilibrium. 
The rotating bar disturbance, with the same strength and 
size as in the perturbative calculation, is again represented 
by the quadrupole term of gravitational potential of the 
ellipsoid. It is turned on adiabatically over four bar ro- 
tation times to avoid sudden transients. The use of the 
quadrupole term only, the first contributing multipole af- 
ter the monopole, ensures that the dark matter halo re- 
mains in approximate equilibrium as the bar perturba- 
tion is applied. The early n-body evolution, shown in 
Figure lb, is similar to the results from the perturbative 
approach. The inner evolution follows the analytic predic- 
tions up to five bar rotation times. At this point, approx- 
imately 30% of the available angular momentum in the 
bar pattern has been transferred to the halo. Subsequent 
evolution is more rapid, presumably due to the non-linear 
response of the near-resonant orbits, although the details 
remain to be investigated. A similar super-linear increase 



in torque was reported by Hernquist and Weinberg (1992). 
The cusp within ^0.5 bar radii is completely removed af- 
ter only 7 bar rotations. At this time the bar pattern has 
lost all of its original angular momentum, the stars in the 
bar have lost ~ 25% of their angular momentum, and the 
disk as a whole has lost only ~8% of its original angular 
momentum, and the disk as a whole has lost only ~8% of 
its original angular momentum. 

3. DISCUSSION 

3.1. Dynamics 

In the previous section we mentioned that one can esti- 
mate the dynamical friction decay time of a rotating bar 
by applying the Chandrasekhar dynamical friction formula 
to a toy model for the bar consisting of two masses in or- 
bit at the same distance from the center of the galaxy 
but at opposite position angles. This physical argument 
is a gross simplification, however. The Chandrasekhar dy- 
namical friction formula is derived by considering the mo- 
mentum transfer of stars gravitationally scattered by a 
traveling body. The problem is much more complicated 
in a halo where individual orbits are quasiperiodic. The 
gravitational wake, i.e. the response of the dark halo to 
the bar, is not as simple as the Chandrasekhar formula 
would lead one to believe. A given orbit may encounter 
the rotating perturbation many times. If the precession 
frequency leads or trails the pattern, the net torque ap- 
plied to this orbit will be zero on average. This appears 
to be in conflict with dynamical friction but really it is 
not. The bar will only receive a net torque if the gravita- 
tional wake either trails or leads the bar pattern and this 
can only occur if some of the orbital actions are changed, 
something that can only occur at or near resonances. At 
exact commensurabilities between the orbital and pattern 
frequencies individual orbits receive kicks breaking sym- 
metry, resulting in a density response that trails the bar. 
To have this symmetry broken requires gradients in the 
phase space density at the positions of these resonances. 
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Fig. 2. — The halo response and bar position during the evolution of an NFW model halo. We subtract the mean density and plot the 
amplitude of the resulting wake where white represents underdense regions and black overdense regions. 



In Figure 2, we plot the actual gravitational wake in- 
duced by a bar in an n-body simulation. It is similar to 
the simulations described above but uses an NFW halo re- 
alized with 10^ equal mass particles and is plotted after 5 
bar pattern rotations. This figure shows the halo density 
distortion in response to the bar. The bar's size and phase 
are shown schematically. The bar pattern leads the wake 
and therefore the bar torques the halo. 

The location in the halo of the peak angular momen- 
tum transfer depends on the halo profile itself in two 
ways. First, the torqued orbits will be near commensu- 
rabilities between the orbital frequencies and the pattern 
speed. These commensurabilities or resonances take the 
form Ir^ir + = mflbar where the three values of f2 

are the radial, azimuthal and pattern frequencies, respec- 
tively, Ir and are integers and m is the azimuthal multi- 
pole index. For the I = m = 2 multipole, we will denote a 
particular resonance by the tuple (Z^, l^). Low-order (high- 
order) resonances have small (large) values of \lr\ or j/^j. 
Second, a net torque requires a differential in phase space 
density on either side of a particular resonance. It is not 
the dark matter density of the inner halo but a large phase 
space gradient that is key to a large bar torque. In fact, 
if the bar were located in an infinitely large homogeneous 
core, no matter how dense, there would be no net torque, 
in contrast to the predictions of the simple toy model. If 
the bar is inside the halo core, then there will be very 
few nearby resonances with a differential in phase space 
density and the torque would be diminished; the domi- 
nant resonant orbits in this case would be at or beyond 
the core radius. Conversely, if the bar is in a dark matter 



cusp, orbits near and inside the bar radius cover a large 
range of frequencies. There are low-order resonances deep 
in the cusp where the phase space gradient is large. The 
torqued cusp orbits move to larger radii, decreasing the 
cusp density and decreasing the overall depth of the po- 
tential well. Both of these effects cause the cusp to expand 
overall. Thus, the formation of a bar will naturally elimi- 
nate an inner cusp. 

Linear theory allows us to explicitly identify the dom- 
inant resonances and the angular momentum transfer to 
halo orbits for models both with and without cores. The 
overall torque is dominated by one or two resonances in 
each case. For the NFW profile, the torque is dominated 
by the resonance (—1,2); this is analogous to an inner 
Lindblad resonance. This resonance does not occur for 
a comparable King (1966) model. The first contributing 
resonance for the King model is (2, —2) near the location 
of the core radius. Higher-order resonances that occur at 
larger radii are well confined to a single characteristic ra- 
dius (close to vertical in Fig. 3. However, the (—1,2) 
resonance is unique. It always exists in a cusp as r ^ 
because fir —>■ 2fl^ as the orbital angular momentum J ap- 
proaches zero. Therefore there is always some value of J 
for which — $7^ + 2570 = 2fli,ar even though fir and fl^ both 
diverge for small r. For this reason, the (—1,2) resonance 
can affect orbits deep within a cusp, dramatically chang- 
ing the inner profile. For a model with a core, the core 
expands somewhat as angular momentum is transferred 
toward the outer halo but otherwise remains qualitatively 
similar to its initial state. 

Since bar instability is ubiquitous in self-gravitating 
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Fic;. 3. — Location of low-order resonances in energy (lower axis) and characteristic radius (upper axis) for the bar in the NFW profile. The 
vertical axis describes the orbital angular momentum J in units of the maximum angular momentum for a given energy, Jmax- The inner 
Lindblad-like resonance extends throughout the inner cusp. This resonance is absent from dark halo models with cores. 



disks, continued accretion and cooling is likely to precipi- 
tate a bar instability in the forming inner gas disk early on. 
Once the rotating bar forms, its body angular momentum 
can be transferred to the halo as described above, flat- 
tening and removing the inner cusp. It is possible for this 
process to occur in stages with multiple bars; at each stage 
the inner core will grow. There will be a transition in the 
magnitude of the torque and a slowing of halo evolution 
when the inner (—1,2) resonance is finally eliminated. 

3.2. Implications 

Because the bar-cusp coupling depends on near- 
resonant dynamics, simulations with very high resolution 
will be required to resolve the dynamics properly. We 

are confident in our n-body results since they agree with 
the exact perturbative results in the linear regime, once 
we have enough particles to reach convergence. To test 
for convergence we ran a suite of simulations with in- 
creasing particle number from 10* through 4 x lO'^ and 
determined that, for our SCF expansion method, we ob- 
tained convergence in the evolution for particle numbers 
> 4x 10^. This class of potential solver (see Glutton-Brock 
1972, 1973; Kalnajs 1976; Fridman and Polyachenko 1984; 
Hernquist and Ostriker 1992) suppresses small scale fluc- 
tuations. Direct summation approaches (e.g. Kawai et al. 
2000, GRAPE), tree codes (Barnes and Hut 1986) and 
grid based codes (Sellwood and Merritt 1994; Pearce and 
Gouchman 1997) have more inherent small scale noise and 
will most likely require even higher particle numbers to 
obtain the same convergence. Note that even if the res- 
onances are not resolved, these simulations will still ex- 
hibit significant torque. Our suite of simulations shows 
that the overall torque increases as the particle number 
decreases^ I These same simulations give us good agree- 

^ If the orbit difTusion becomes large enough, it is plausible that the torque will decrease to zero, however. 



ment with Chandrasekhar's formula. Such agreement is 
not a good indication that one is observing the correct 
dynamics. Gonversely, the Ghandrasekhar formula works 
well in simulations with low to moderate particle number 
because the resonances are obliterated by artificial diffu- 
sion and, therefore, are well represented by simple scatter- 
ing. In short, it is difficult to see resonant effects in n-body 
simulations because the diffusion rate is high for moderate 
numbers of particles. Astronomical sources of noise such 
as orbiting substructure, decaying spiral waves, lopsided- 
ness, etc. do not produce enough small scale noise to affect 
this resonant evolution but instead produce large scale de- 
viations from equilibrium that will not drive significant 
relaxation in the inner halo within several orbital times 
(see Weinberg 2001b, a, for estimates of these timescales). 

The evolution of bars within dark halos has recently 
been studied by Debattista and Sellwood (1998, 2000). 
They simulated bars in dark halos by constructing a strong 
bar from a Q-unstable disk and then following its evolu- 
tion. They found that the bars rapidly transferred an- 
gular momentum to the massive (non-rotating) halos as 
predicted by Weinberg (1985). On the other hand, a few 
direct and a wider variety of indirect inferences from ob- 
servations indicate that galaxies today harbor rapidly ro- 
tating bars. Hence, Debattista & Sellwood concluded that 
dark matter halos must have low density (or large cores) to 
be consistent with observations and could not have cusps. 
These results are consistent with our scenario. As the disk 
matures and becomes stellar rather than gas dominated, 
a normal stellar bar may form through secular growth or 
instability. The first generation of bar evolution would 
have eliminated the inner halo torque by removing the 
cusp and, hence, be consistent with the Debattista and 
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Sellwood (2000) arguments. 

The bar torque will affect the kinematics of the bulge. 
Kormendy (1982) finds that bulges of SB galaxies have suf- 
ficient spin to be rotationally flattened. However, triaxial 
galaxy bulges appear to be kincmatically midway between 
an isothermal and rotationally-supported disk component. 
The bulge length scale is smaller the classic bar length and 
much smaller than the length of a strong primordial bar. 
Because the bulge profile is not supported by relatively 
low-energy eccentric orbits, the bulge will not be subject 
to the strong density evolution predicted here. 

Most of the evidence against central dark matter cusps 
in galaxies concerns dwarf galaxies, particularly those with 
low surface brightness (Cote et al. 2000; de Blok et al. 2001; 
Blais-Ouellette et al. 2001) since these systems offer more 
precise rotation curve decompositions. One might not ex- 
pect strong bars to form in such systems and, hence, for 
the mechanism proposed here not to have much relevance. 
However, the same analysis used to indicate the lack of a 
central dark matter cusp also shows that these low sur- 
face brightness galaxies are deficient in baryons by three 
or four times compared to normal galaxies (van den Bosch 
and Swaters 2001). If a strong bar forms in a gas rich disk 
of a dwarf galaxy, much of the gas will lose substantial 
amounts of angular momentum and be driven towards the 
center (Roberts et al. 1979). This dense cold gas would 
then be expected to undergo a strong starburst and, ow- 
ing to the shallow potential well of the dwarf galaxy, much 
of the gas would be expelled as a wind (Dekel and Silk 
1986). The remaining galaxy would be one of low surface 
brightness possessing a core in its dark matter distribu- 
tion. Our bar mechanism, therefore, not only provides a 
natural explanation for the existence of dark matter cores 
but for the existence of low surface brightness dwarfs as 
well; those dwarf galaxies that had the strongest bars will 
have the largest cores and lowest surface brightness. There 
is also some evidence that larger, high surface brightness 
galaxies also have cores in their central dark matter dis- 
tributions, although such determinations are much more 
difficult (Debattista and Sellwood 2000; Binney and Evans 
2001). 

Even if our proposed bar mechanism does not explain 
the lack of observed central dark matter cusps in all cases. 



it will still have a profound effect on the structure and evo- 
lution of almost all galaxies. Hence, both to understand 
galaxy formation and evolution and to make predictions 
from theory it is necessary to resolve these dynamical pro- 
cesses. For example: it changes the predicted dark matter 
densities at the solar radius and would change the pre- 
dicted dark matter detection rates (see Stiff et al. 2001, 
and references therein), it could change the rotation ve- 
locities used in TuUy-Fisher predictions and perhaps im- 
prove the agreement of theory with observations (Navarro 
and Steinmetz 2000); it could change disk scale lengths 
(Steinmetz and Navarro 1999); and it could make satellite 
galaxies easier to tidally disrupt by decreasing their central 
densities (Moore et al. 1999). Some of these consequences 
are qualitatively discussed in Binney et al. (2001) but it is 
hard to address many of these issues quantitatively with 
the simplified simulations presented here. In the future, 
we plan to perform simulations with self consistent disks 
and bars, eventually including gas dynamics, to investigate 
these issues. 

The lack of an observed central dark matter cusp in 
galaxies is a consequence of simple dynamical evolution 
and does not require a fundamental change to the nature 
of dark matter or galaxy formation. The difference in evo- 
lutionary end states may be the result of strong star for- 
mation feedback evacuating the more weakly bound cen- 
tral potentials, lack of strong accretion events and mergers 
after the primordial bar has disappeared, or a combina- 
tion of the two. Indeed, such stochastic processes natu- 
rally predict the inferred dispersion in present day profiles. 
However, to correctly resolve these important dynamical 
processes in ab initio calculations of galaxy formation re- 
mains a daunting task, requiring at least 4,000,000 halo 
particles using our SCF code, and requiring many times 
more particles when using noisier tree, direct summation, 
or grid based techniques — the usual methods employed in 
such calculations. 
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